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Abstract 

Biological oscillators coordinate individual cellular components so that they function coherently 
and collectively. They are typically composed of multiple feedback loops, and period mismatch is 
unavoidable in biological implementations. We investigated the advantageous effect of this period 
mismatch in terms of a synchronization response to external stimuli. Specifically, we considered 
two fundamental models of genetic circuits: smooth- and relaxation oscillators. Using phase reduc- 
tion and Floquet multipliers, we numerically analyzed their entrainability under different coupling 
strengths and period ratios. We found that a period mismatch induces better entrainment in both 
types of oscillator; the enhancement occurs in the vicinity of the bifurcation on their limit cycles. 
In the smooth oscillator, the optimal period ratio for the enhancement coincides with the experi- 
mentally observed ratio, which suggests biological exploitation of the period mismatch. Although 
the origin of multiple feedback loops is often explained as a passive mechanism to ensure robustness 
against perturbation, we study the active benefits of the period mismatch, which include increasing 
the efficiency of the genetic oscillators. Our findings show a qualitatively different perspective for 
both the inherent advantages of multiple loops and their essentiality. 
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1. INTRODUCTION 



Biochemical oscillators exist primarily at the genetic level, such as in the interplay among 
mRNAs and proteins At least one negative feedback loop (NFL) is essential for limit- 
cycle oscillation, and in the minimum configuration, one NFL with time-delay is sufficient 
to generate oscillation. However, many realistic oscillators possess extra components, such 
as positive feedback loops (PFL) and more than one NFLs, that are apparently redundant. 
Since more regulatory components consume more resources (e.g., amino acids for synthesis 
and adenosine triphosphate for operation), this would seem energetically disadvantageous, 
and thus not preferable through evolution. In addition, the coexistence of positive and 
negative loops may make the resulting behavior incoherent. Still, the circadian clock of even 
the simplest life form (e.g., the marine dinofiagellate Gonyaulax polyedra)is composed of 
multiple oscillators whose rhythms exhibit different periodic patterns 0, 3]. Usually this 
observation is explained by the notion of biological robustness, in which there is a backup 
mechanism or a resistance to perturbation. However, such observations do not suggest the 
essentiality of multiple loops: why does not even a single species possess the minimum 
configuration? In this study, we numerically analyze the inherent advantages of multiple 
loops from the perspective of entrainment and the adaptation of the system to external 
input 

Biologists have investigated coupled oscillators at the molecular level in many different 
organisms. Bell-Pedersen et al. reviewed circadian clocks in several species from different 

n 

clades [5|. They suggested that multiple loops comprise pacemaker and slave oscillators 
in unicellular organisms, whereas in mammals and avians, a centralized pacemaker entrains 
downstream systems. Circadian systems have been also investigated theoretically. Wagner et 
al. studied the stability of an oscillator (called the Goodwin model [g]) against perturbation 
of the kinetic rate, and found that interlinked loops are more robust (by nearly one order 

n 

of magnitude) than a non- linked counterpart [7]. Other benefits of multiple feedback loops 
are found in Refs , and their intuitive advantages are summarized in the review by 
Hastings Q. 

When multiple loops are involved, synchronization becomes a crucial issue. In biological 
implementations, no loops can share an identical period in a strict sense. At first glance, 
such disorder seems to degrade the system performance. However, the results in nonlinear 
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physics indicate that a certain amount of disorder actual 



mechanisms similar to those in stochastic resonance 
noise-induced synchronization 
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y enhances performance through 



Ill4l8l|. Brownian transport 19l-|23|. 



and disorder- induced resonance j27| . 



Our focus is the property called entrainability, or the ability of oscillators to be synchro- 
nized with an external periodic signal. It is directly connected to the adaptability of the 
species and therefore its survival l28| . Not only does a Drosophila with a circadian clock 
survive better than one without ^91, but also an Arabidopsis with better entrainment with 

n 

the environment has increased photosynthesis and consequently grows faster |30j. Even 
cyanobacteria with better entrainability have an advantage over nonsynchronizable ones, 
although this advantage disappears in constant environments 31 |. 

At the molecular level, entrainability has been attributed to several different mechanisms. 
Gonze et al. studied a population of suprachiasmatic nuclei (SCN) neurons and showed that 
efficient global synchronization is achieved via damping 32|. Liu et al. experimentally 
confirmed that robust oscillatory behavior can be achieved by tissue-level communication 
between SCN cells, even if critical clock genes have been knocked out 33|. Tsai et al. 
reported that the presence of PFLs allowed for easier tuning of the period without changing 



the amplitude 



34| . Zhou et al. showed that positive feedback loops assist noise-induced 



synchronization through their sensitivity to external signals 35 1. 

Still, a fundamental question has not been answered: what is the advantage of harboring 
multiple loops with mismatched periods? We investigate the effects of mismatched periods 
on limit-cycle oscillations from the viewpoint of entrainment. Specifically, we consider the 
minimal fundamental model: an oscillator made of two loops, each of a different period, that 
are connected. By using phase reduction and Floquet multipliers, we analyzed the response 
of this system to external periodic stimuli as a function of the coupling strength and the 
period ratio. Our main finding is the large effect of mismatched periods on entrainability: 
that is, better entrainment is achieved when the two oscillatory components have different 
natural periods. This enhancement effect is observed in two types of coupled oscillators: 
smooth and relaxation oscillators (Figured]). The enhanced entrainability seen in the present 
study is similar to that reported by Komin et al. [36]. However, our mechanism is different 
from theirs: the enhancement in our case is caused by the oscillators in the vicinity of the 
bifurcation on the limit cycle whereas their enhancement results from damped oscillation 
(oscillator death). 
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2. METHODS 



We studied the effects of mismatched periods on entrainment by using smooth and relax- 
ation oscillators. Before explaining the details of these models, we will first explain period 
scaling and a coupling scheme. 

2.1. The base model 

An A^-dimensional differential equation for a single oscillator is represented by 

l^/W, (1) 

where v and f{v) are A^-dimensional column vectors defined by 

v = {x,y,z,--- y, f{v) = f2{v), ■■■ , fN{v)y, 

where T denotes a transposed vector (matrix). Equation [T] yields autonomous oscillation (a 
limit-cycle oscillation) whose period is T. The period can be tuned without changing the 
amplitude of the oscillation by introducing a scaling parameter r 

dv _ f{v) 

dt~ T ' 

where the period of equation |2] is given by rT. Other than their periods, solutions to 



equation [2] with different r values will have the same properties [32|, l36|. We used equation [2] 
to study the effects of mismatched periods between two oscillatory components: 

dvi f{vi) 



dt Ti 

dV2 f{V2] 



Ci{v,,V2), (3) 

+C2{V2,Vi), (4) 

at T2 

where the subscript i denotes the ith component and Ci{vi, Vj) represents a coupling term. 
To investigate the effects of mismatched periods, we set 

n = 1, T2 = R, (5) 

where R represents the ratio of the periods of the two oscillators. 

2.2. Genetic oscillator models 

We show two specific models of f{v), each for the smooth and relaxation oscillators. 
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FIG. 1: (A) Diagram of a smooth oscillator, (B) simplified diagram of a single smooth oscillator, 
(C) two coupled smooth oscillators, (D) diagram of a relaxation oscillator, (E) simplified diagram 
of a single relaxation oscillator, and (F) two coupled relaxation oscillators. The symbols — >■ and H 
represent positive and negative regulations, respectively. 
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FIG. 2: Time courses of (A) the smooth oscillator (equation [6j) and (B) the relaxation oscillator 
(equation [8]) . Solid, dotted, and dashed lines denote x, y, and z, respectively. The period and the 
amplitude of y are normalized to 1 in both models. 



2.2.1. Smooth oscillator 

A smooth oscillator is the most basic structure for limit-cycle oscillations and many 
biochemical oscillators (e.g., the Goodwin model ^ and the repressilator js?!) belong to 
this class. We adopted the smooth oscillator presented in Novak and Tyson [l] (Figure [1]). 
A diagram and the simplified structure of the smooth oscillator are shown in Figures [TJ4. 



and B, where the transcription of X (mRNA) is inhibited by Z (protein). Protein Y plays 
a role in the time-delay of the loop X ^ Y ^ Z -\ X , where X Y and X -\ Y 
denote activation and repression, respectively. In the presence of sufficient time delay in 
the NFL, the regulatory network exhibits limit-cycle oscillations. The dynamics of the 
smooth oscillator are described by the following three-dimensional differential equation on 
V = {x,y,z)~^, where x, y, and z denote the concentrations of mRNA X, protein Y, and 
protein Z, respectively: 



hgyX 



y 



dy- 



\ 



(6) 



Km + y 

Because the amplitude and the period of the limit cycle can be tuned separately by proper 
scaling, we employed the parameters = 5.44, kdx = 8.99, kgz = 4.49, kdz = 4.49, 
ksy = 18.0, Vdy = 2.72, Km = 0.00303, Kd = 0.303, and p = 4.0 to achieve unit period and 
unit amplitude for y (we define the amplitude as the difference between the maximum and 
the minimum of the oscillation. See Figure [2JA.) . These parameter settings are essentially 
identical to those presented in the Supplement of Novak and Tyson [l| , except for the scale 
transform (V^, = 0.2, kd^ = 0.1, ksz = 0.05, kdz = 0.05, ksy = 0.2, Vdy = 0.1, = 0.01, 
Kd = 1.0, and p = 4.0 in the original model). 

We analyzed two coupled smooth oscillators, as shown in Figure [Tp, where the subscript 
i of Xi denotes the ith oscillator component. We assumed that the protein translated from 
Xi positively regulates the transcription of X2. The coupling term in equations [3] and H] can 
be written in the following linear relation: 

\ 

((z,j) = (l,2)or (2,1)), (7) 



Ci{vi,Vj 



\ 









where e denotes the coupling strength. 



2.2.2. Relaxation oscillator 



The existence of PFLs makes oscillators excitatory, and such oscillators are known as 

n 

relaxation oscillators [JJ. Relaxation oscillators have cusps in their nuUcline curves, and one 
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example of this class is the FitzHugh-Nagumo oscillator 38|, |39 1 . 



We again employ a relaxation-type genetic oscillator presented by Novak and Tyson 
Figures [Tp and E show specific and simplified models, respectively, of a relaxation oscillator 
in which transcription of X (mRNA) is inhibited by Y (protein), which serves as a NFL. 
Following this, Y (mRNA) is degraded by Z (enzyme) and also forms an enzyme-substrate 
complex ZY to inhibit the functionality of Z {ZY + Y ^ ZYY). The interaction Y -\ Z -\Y 
functions as the PFL in this network. In many cases, it is assumed that the concentration of 
enzyme Z equilibrates (i = 0), yielding the following two-dimensional differential equation 
in terms of i; = {x,y)^, where x and y denote the concentrations of mRNA X and protein 
Y, respectively: 



fiv) 



(8) 



As in the smooth oscillator, we scaled the parameters so that the period and the amplitude 
y both became 1 (Figure OB). The parameters, Vgx = 0.815, = 0.0268, k^x = 3.04, 
ksy = 60.8, kdy = 3.04, Vdy = 16.3, Ki = 0.134, Kd = 0.268, and p = 4.0, are essentially 
equivalent to those presented in the Supplement of Novak and Tyson [L] , except for the scale 
transform {V^x = 0.05, = 0.1, kdx = 0.05, k^y = 1.0, kdy = 0.05, Vdy = 1.0, = 0.5, 
Kd = 1.0, and p = 4.0 in the original model). 

We analyzed two coupled relaxation oscillators, as shown in Figure HF. The coupling term 
C,' was written ClS db linear relation: 

C.{v„v,) = I ((,,^-) = (1,2) or (2,1)), (9) 

where e denotes the coupling strength. 



2.3. Analytical approaches 

For an analytical analysis, we used phase reduction and Floquet multipliers. Even a 
system of only two coupled oscillators exhibits very complicated behaviors, such as n : m 
synchronization (n, m; positive integers) and chaos, depending on the model parameters (see 
section |3]). Because detailed bifurcation analysis of the coupled oscillator is outside the scope 
of this paper, we assumed 1 : 1 synchronization between two oscillatory components (tight 
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coupling). Under this assumption, the coupled oscillator can be viewed as one oscillator, 
whose dynamics is represented by 

du , . , , 

* = (10) 



with 

(11) 



u = \ M , g{u) 



V2 



y — + C2{V2,Vi) J 



T2 

Incorporating an input signal in equation [TOl we have 

^=g{u) + xliujt), (12) 

where I{ut) denotes an input signal with angular frequency u and x is a scalar representing 
the signal strength. We will use equation [12] in the calculations for entrainability and Floquet 
multipliers. 

2.3.1. Phase reduction 

A phase-reduction approach uses a phase variable to reduce a multidimensional system 



into a one-dimensional differential equation 40|, |4l| . In the absence of external perturbations 



(equation [TOj) , a limit-cycle oscillation can be described in terms of a phase variable G 
[0,27r]: 

where Vt is the natural angular frequency of the limit cycle {Vt = 2tt/T, where T is the 
period). Although the phase in equation [13] is only defined on the stable limit-cycle 
trajectory, its definition can be expanded to the entire u space, where the equiphase surface 
is called the isochron X(0). An example of an isochron X(0) defined on a hypothetical limit- 
cycle oscillator, shown by dashed lines, in Figure [31^ where the isochron is drawn at intervals 
of 7r/6. 

In the presence of external signals (equation [T2|) . the time evolution of (f) is described by 

^ = fi + xC/(0)-JM, (14) 

where a dot (■) represents an inner product and C/(0) is the phase response curve (PRC) 
defined by 

t/(0') = V0(n)|_^(^,). (15) 
8 
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FIG. 3: (A) Trajectory of a hypothetical hmit-cycle oscillator (solid line) and its isochron 
(dashed lines), where the isochron is drawn at intervals of 7r/6. (B) Illustration oi = + 
^Q(^) = 0, where the solid and dashed line denote X0(V') aiid —Ail, respectively, for typical 
cases. A stable solution (empty circle) exists inside x0max and x0min, whose length is denoted by 
xW . (C) Illustration of the relation between the Arnold tongue (colored region), inside which the 
oscillator can synchronize to an external signal (vertical and horizontal axes represent the signal 
angular frequency uj and the signal strength x, respectively). Equation [19] (dashed line) is a linear 
approximation of this region for sufficiently small x- The width of the region can be approximated 
by x(0max - ©min) around X = 0. 

Here Wo(0') represents a point on limit-cycle oscillation at the phase 0'. It is assumed that 
the perturbed trajectory is in the vicinity of the unperturbed limit-cycle trajectory (i.e., x 
is sufficiently small). The introduction of a slow variable = (p — ut m equation \\M yields 



^ = Ail + xUi^J + ^i) ■ I{^t), 



(16) 



where Afi = Vt — u. Because tp is a, slow variable {ip is very small), when x is sufficiently 
small, the inner product term can be approximated by its average (separating the timescales 
of tp and U{%1) + out) ■ I{ujt)): 



^ ^ AQ + xe{^), 



with 



2tt 



Q(^) = — / deui^ + e)-iie), 



(17) 



(18) 



(the integration of equation [18] is evaluated numerically in practical calculations). If stable 
solutions exist, this oscillator synchronizes to external signals. Q{iIj) is a 27r-periodic func- 



tion, and its shape is typically given by sinusoidal- like functions, as shown in Figure |3j3. The 
fixed points are the intersection points of and —AQ, and only the point with ^ < 
is a stable solution (an empty circle in Figure [3B)- Therefore, the oscillator is synchronized 
to external signals if the following relation holds: 

+ n<cu< xQm.. + ^, (19) 

with 

0mm = mine(V'), Gmax = max e 

If) if) 

The extent of synchronization in terms of x (signal strength) and u (signal angular fre- 
quency) is often described by a domain known as the Arnold tongue or the synchronization 
region 4^, inside of which the oscillator can synchronize to an external periodic signal 
(the colored domain in Figure [3]C). Although the border of the Arnold tongue is generally 
nonlinear as a function of the upper and lower borders can be approximated linearly by 
equation [T9] when x is sufficiently small (these are shown with the dashed line in Figure [3]C). 
Thus the range in which u can be synchronized with an external signal can be quantified 
by x(0max — ©mill), which represents the width of the Arnold tongue for sufficiently small 
X- We define the entrainability as a x-independent variable (see Figure [3l3): 

w = e^ax - e^i,. (20) 

Writing equation [19] in terms of the period of the external signal, we have 

<Ts< -7^ (21) 



where is the period of the external signal [Tg = lnjiS) (note that equations [T9] and [2T] 
only hold for sufficiently small x). The PRC [equation [15] can be calculated by the following 



differential equation |4ll, |43| : 



^ = -{Dg(no(t))rC/(t), (22) 

with the constraint 

t/(t)-g(no(t)) = fi, (23) 

where u^it) is a stable periodic solution of equation [TD] iu^it + T) = ui^(t), where T is 
the period), Dg{uo{t)) is a Jacobian matrix of g{u) around UQ{t), and we abbreviated 
U(t) = U{(j)(t)) {(pit) is the phase as a function of t according to equation [T3|) . Because 
equation [22] is unstable, we numerically solved equation [22] backward in time. 
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2.3.2. Floquet multiplier 



In the presence of external signals, the solution u{t) deviates from a stable orbit no, 
where the deviation r) = u — uq obeys the differential equation 

^ = Dg{u,{t)Ut). (24) 

Here we assume that the deviation is sufficiently small. Because Dq^Uq) is a periodic 
function with period T, equation |2l] has the Floquet solutions, generally represented by 

N 

r]{t) = Y^ a exp{^it)pi{t), (25) 
1=1 

where Pi{t) are functions with period T, Cj is a coefficient determined by the initial conditions, 
and /ij are the Floquet exponents. The Floquet multipliers pi are related to the exponents 
via 

pi = exp(/iiT). (26) 

Note that the Floquet exponents are a special case of the Lyapunov exponents for periodic 
systems (e.g., limit cycles) [4^. The Floquet multipliers can be calculated numerically from 



a matrix differential equation in terms of the matrix Q{t), after Klausmeier 44j |: 

^ = Dg{uo{tmt). (27) 

Equation [27] is integrated numerically from t = to t = T with an initial condition that 
Q{t = 0) is the identity matrix. The Floquet multipliers are obtained as the eigenvalues of 
Q{t = T), and hence A^- dimensional limit-cycle models yield N Floquet multipliers. 



When the Floquet multipliers pi cross the unit circle \p\ = 1 on the complex 
limit cycle bifurcates. Bifurcations of codimension 1 can occur in three patterns 



jlane, the 



45|-|48|: 



• Saddle-node bifurcation: A single real Floquet multiplier crosses p = 1 (this includes 
transcritical or pitchfork bifurcations as special cases); 

• Hopf bifurcation: Complex-conjugate Floquet multipliers cross \p\ = 1 (p 7^ ±1); 

• Period-doubling bifurcation: A single real Floquet multiplier crosses p = —1. 

In autonomous systems, one of the Floquet multipliers is = 1, which corresponds to a pure 
phase shift. Wiesenfeld and McNamara showed that Floquet multipliers whose absolute 
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values are close to 1 are responsible for signal amplification j47|, |48|. We calculated the 
Floquet multipliers up to the second-largest absolute value (except for the inherent multiplier 
p = 1), which are denoted as pi (the largest) and p2 (the second largest) in our models. 



3. RESULTS 

We performed a bifurcation analysis and examined the relationship between entrainability 
and the Floquet multiphers. We set ti = 1 and T2 = R in equations E] and HI where R 
represents the ratio of the periods of the two oscillators. When R ^ 1, the two oscillatory 
components have different natural periods (without coupling) and hence the periods are 
mismatched. As the controllable parameters for the model calculations, we adopted R and 
e, where e denotes the coupling strength defined in equations [7] and [91 

3.1. Bifurcation diagram 

Bifurcations of the limit cycle are classified into bifurcations in equilibria and those on 
the limit cycle. Here, bifurcation on the limit cycle was studied by numerically solving 
differential equations El and El 

Specifically, we described the bifurcation diagrams as a function of R (period ratio) and 
e (coupling strength) and calculated the local maxima of yi. If the oscillator showed regular 
limit-cycle oscillations, the maxima of yi was represented by a single point. If the limit cycle 
underwent bifurcation, the maximal points were described by more than two points. We 
iterated by changing R from -Rmin to -Rmax (-Rmin < Rmax) and then adopted the last values 
of the preceding R values as the initial values for the next R value. The incremental step 
for R was set at AR = (-Rmax — -Rmin)/1000, where -Rmin and -Rmax are indicated as the lower 
and upper limits on Figures HI and [61 At each parameter value, we discarded as transient 
phases those during t G [0,200], and we recorded the local maxima of the trajectories in 
during the intervals t G [200,300]. We carried out these same procedures on the coupling 
parameter e (emin, Cmax, with the increment Ae = (emax — emm)/1000). 

In Figure [H the bifurcation diagram of yi as a function of -R in the smooth oscillator is 
shown for (A) e = 1.5 and (B) e = 2.5. Although the maxima of yi slightly depend on -R, 
there is no bifurcation. In Figure [5l the bifurcation diagram as a function of e is shown for 
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FIG. 4: Bifurcation diagram (local maxima of yi) of the smooth oscillator as a function of i?, where 
(A) e = 1.5, and (B) e = 2.5. 

(A) R = 1.4 and (B) R = 2.5. In both cases, bifurcations occurred a.t e 1.3 for R = 1.4 
and e ~ 0.9 for R = 2.5, and chaotic behavior was observed when the coupling strength was 
small. 

Figure [6] shows the bifurcation diagram as a function of R for the relaxation oscillator 
for (A) e = 2.0 and (B) e = 2.5. In the former case, a period-doubling bifurcation occurred 
in the region 1.7 < R < 2.7. In the latter, a cusp appeared around R = 1.9, in the vicinity 
of a saddle-node bifurcation (see the section on entrainability and Floquet multipliers). In 
Figure[7l the bifurcation diagram as a function of e is shown for (A) R = 1.3 and (B) R = 1.9. 
A complex behavior (chaos) appeared in a small region of e. Especially in Figure [7]B, we 
found a branch of a maxima around 1.7 < e < 2.3, indicating the occurrence of the period- 
doubling bifurcation. 

In comparison, bifurcations are less likely in a smooth oscillator than in a relaxation 
oscillator; bifurcations are only observed in the case of very weak coupling in smooth oscilla- 
tors (Figure [5]). In contrast, the relaxation oscillator exhibited period-doubling and chaotic 
oscillations in a wider range of coupling strengths (Figures [6] and [7]). 

Next, we examine the time course of the coupled oscillators with mismatched periods 
(Figures [8] and [H respectively, for the smooth and relaxation oscillators). Without mis- 
matched periods, the temporal trajectory of the oscillators would be identical to that of the 
case of a single oscillator (Figure [2]). Figure [8] shows the case for coupled smooth oscillators 
with (A) e = 1.5, R = 1.5, and (B) e = 2.5, R = 1.5; Figure [9] shows the case for coupled 
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FIG. 5: Bifurcation diagram (local maxima of yi) of the smooth oscillator as a function of e, where 
(A) R = 1.4, and (B) R = 2.5. 




FIG. 6: Bifurcation diagram (local maxima of yi) of the relaxation oscillator as a function of R, 
where (A) e = 2.0, and (B) e = 2.5. 

relaxation oscillators with (A) e = 2.0, R = 2.8, and (B) e = 2.5, R = 1.9. In all the cases 
shown in Figures |8] and HI we see the 1 : 1 synchronization between the two oscillators. 
Without period mismatch, the amplitude of oscillation y is 1 (see Methods). We see from 
Figures |8] and [9] that the amplitudes of yi and 1/2 depend on the period ratio R, where the 
amplitude of yi is smaller than 1 except in Figure [9|3. From Figures H] and El which describe 
maxima of yi as a function R, we see that yi is larger than 1 in both oscillator models when 
the period mismatch is smaller, and the maxima of yi start to decrease with increasing R. 
The decrease of the amplitude is larger for the smooth oscillator than for the relaxation 
oscillator. 
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FIG. 7: Bifurcation diagram (local maxima of yi) of the relaxation oscillator as a function of e, 
where (A) R = 1.3, and (B) R = 1.9. 
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FIG. 8: Time course of the coupled smooth oscillator with (A) e = 1.5 and R = 1.5, and (B) 
e = 2.5 and R = 1.5. 
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FIG. 9: Time course of the coupled relaxation oscillator, with (A) e = 2.0 and R = 2.8, and (B) 
e = 2.5 and R = 1.9. 
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3.2. Entrainability and Floquet multiplier 



In using phase reduction and Floquet multipliers to calculate how the entrainability 
W (equation depends on R (equation E]) and e (equations [7] and [H]), we employed the 
following sinusoidal input signal: 



I{ujt) 



^ sm(ujt) ^ 





(28) 



V '■ J 

where the input signal only affects the concentration of mRNA in the first oscillator (xi in 
equation [T0|) . We assumed 1 : 1 synchronization between the two oscillators (the regions 
where the maximum of yi is represented by a single point in Figures HHTj) so that the coupled 
oscillator could be viewed as a single case. Although we calculated the entrainability and 
the Floquet multipliers inside the 1 : 1 synchronized regions, the entrainability very close to 
the bifurcation points is not shown for some parameters, because equation [22] did not exhibit 
stable periodic solutions even after remaining there for a long time (e ~ 1.3 in Figure [TT1\. 
R ^ 2.7 in Figure fT2]A and e 2.3 in Figure fT3D). In most cases we calculated the Floquet 
multipliers up to the second most dominant ones. When the multipliers were given by 
complex-conjugate values, we checked the three multipliers pi, pi, and p2, where p is the 
complex conjugate of p. When two real multipliers (pi and p[) eventually emerged due to 
annihilation of the complex-conjugate multipliers (pi and pi) by changing parameters, we 
checked the three multipliers pi, p[, and p2. 

By calculating the entrainability W and the Floquet multipliers pi, we obtained stable 
limit-cycle trajectories in a way similar to what we did in the bifurcation analysis. We 
iterated over R and adopted the values of the preceding R values as the initial values for the 
next R values. For each parameter value, we calculated the stable limit-cycle trajectories 
for one period, and these trajectories were used for the calculation of the PRC functions 
(Equation [T5l). These procedures were also carried out for the coupling parameter e. The 
Floquet multipliers also require stable limit-cycle trajectories, and we calculated them in 
the same way as described above. 
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FIG. 10: Entrainability W, period T, and Floquet multipliers pj, as a function of R in the smooth 
oscillator. (A-C) R dependence of (A) W and T, (B) pi, and (C) p2, for e = 1.5. In (A), the 
solid and dashed lines denote the entrainability and the period, respectively. In B and C, the solid, 
dashed, and dot-dashed lines denote the real part, the imaginary part, and absolute value of the 
Floquet multipliers, respectively. 




FIG. 11: Entrainability W, period T, and Floquet multipliers pi, as a function of e in the smooth 
oscillator. (A-C) e dependence of (A) W and T, (B) pi, and (C) p2, for R = 1.4. The meanings 
of the lines are the same as in Figure [TOj 

3.2.1. Smooth oscillator 

We first investigated the smooth oscillator as a function of R (Figure [T0|) . Specifically, 
Figure [TOl A is a dual-axis plot of the entrainability W (solid line; left axis) and the period 
T (dashed line; right axis) for e = 1.5. We can uniquely define the period T of the coupled 
oscillator because we are considering a 1 : 1 synchronization between the two oscillators. In 
Figure fTOl\. the period T is not strongly affected by period mismatch and is around 1 in the 
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whole R region, although the entrainability W depends on R. The entrainability is W ~ 2 
when the two oscillators have the same period [R = 1), whereas ~ 12 for R > 1.4. The 
first and the second dominant Floquet multipliers for e = 1.5 are shown in Figure [TUB and 
C, where solid, dashed, and dot-dashed lines respectively denote Re(pj), Im(pj), and \pi\ 
{pf. Floquet multiplier defined by equation f2E\i . In Figure [TUB, symmetric imaginary curves 
represent complex-conjugate multipliers, and a branch of the real part around R ~ 2.2 is due 
to annihilation of the complex-conjugate multipliers (two real multipliers emerge in place of 
the complex-conjugate multipliers). We see that the pi in Figure [TUB are imaginary values, 
and their absolute values are close to 1 for the i? > 1.4 region. The coupled oscillator is 
close to the Hopf bifurcation points. The second dominant Floquet multiplier p2 is described 
in Figure [TU1C. and we see that its magnitude is comparable to that of pi. The coupled 
oscillator is near the saddle-node bifurcation points (the second Floquet multiplier p2 is a 
pure real value). We also carried out the same calculations with e = 2.5 and found enhanced 
entrainability ~ 16 in the presence of a period mismatch {R > 1.5), as in Figure [TOlC 
(data not shown). According to the analysis with Floquet multipliers, this enhancement can 
also be accounted for by the bifurcation. 

Figure [TT] shows the entrainability and Floquet multipliers as a function of e in the smooth 
oscillator while keeping R constant. Figure [TTK is a dual-axis plot of the entrainability W 
(left axis) and the period T (right axis) for R = 1.4, where the period T does not strongly 
depend on the coupling strength e. In contrast, the entrainability W achieves a maximum, 

~ 15 at e ~ 2.5, which is more than 7 times larger than the no-mismatch case {W ~ 2). 
Figure [TTB shows the first dominant Floquet multiplier pi as a function of e, where the 
notation is the same as in Figure [TUl pi is a complex multiplier for e < 5.5, and its 
absolute value (dot-dashed line) approaches 1 as e decreases, showing that the oscillator 
is in the vicinity of the Hopf bifurcation. We see consistent results in Figure [5] where 
the oscillator is quasiperiodic (chaotic) in the e < 1.3 region. At e ~ 5.5, the complex- 
conjugate multipliers annihilate, and two real multipliers emerge instead. Figure HVC shows 
the second dominant Floquet multiplier p2 as a pure positive real number, responsible for 
the saddle-node bifurcation points (this multiplier is close to pj = 1). We see that p2 achieves 
its maximal value around e ~ 2.5, where the entrainability W also exhibits a maximum. 
Because |pi| in Figure [TTB monotonically decreases as a function of e for e < 5.5, pi does 
not explain the qualitative behavior of W, whereas p2 does. This result shows that being 
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FIG. 12: Entrainability W, period T, and Floquet multipliers function of R for the 

relaxation oscillator. (A-C) R dependence of (A) W and T, (B) pi, and (C) p2 for e = 2.0. (D-F) 
R dependence of (D) W and T, (E) pi, and (F) p2 for e = 2.5. The meanings of lines are the same 
as in Figure [lOl 

close to the saddle- node bifurcation is important for achieving entrainment. We investigated 
the entrainability with a different parameter for the period mismatch R = 2.5, and observed 
that ^ 26 was maximized at an intermediate coupling strength e ~ 5 (data not shown). 
The Floquet multiplier analysis also indicated that enhancement in the R = 2.5 case is 
caused by the oscillators being in the vicinity of the saddle-node bifurcation on the limit 
cycle. 

3.2.2. Relaxation oscillator 

We next apply the same analysis to the relaxation oscillator: we calculated the entrain- 
ability ly, the period T, and the Floquet multipliers pi, as functions of R (Figure [12]) and 
e (Figure [13]). Figure [T2K shows a dual-axis plot of the entrainability W and the period 
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FIG. 13: Entrainability W, period T, and Floquet multipliers pi, as a function of e for the relaxation 
oscillator. (A-C) e dependence of (A) W and T, (B) pi, and (C) p2, for R = 1.3. (D-F) e 
dependence of (D) W and T, (E) pi, and (F) for R = 1.9. The meanings of the lines are the 
same as in Figure [TOl 

T, for e = 2.0 (see also Figure fTOll . The coupled oscillator is not 1 : 1 synchronized in the 
region 1.7 < R < 2.7, and these unsynchronized regions were not plotted. Since the entrain- 
ability W grows rapidly in the vicinity of the bifurcation points {R ~ 1.7) in Figure [T2j\, 
we calculated the first pi and second p2 dominant Floquet multipliers in Figure [T2B and C, 
where pi is a complex-conjugate multiplier around R = 1.7 and is a pure real multiplier very 
near the bifurcation points. We see that pi of e = 2.0 exhibits a sudden decrease towards 
— 1 near R ~ 1.7 and R ^ 2.7, indicating that the coupled oscillator is close to the period- 
doubling bifurcation points. The bifurcation diagram of Figure [6] shows a branch, i.e., the 
period-doubling bifurcation, pi is a real multiplier for R < 3.3, and its value approaches 
— 1 as R decreases toward R = 2.7 in the vicinity of the period-doubling bifurcation. At 
R ~ 3.3, complex-conjugate multipliers emerge instead of the annihilation of two real mul- 
tipliers. Figure \TTC describes the second dominant Floquet multiplier p2 of e = 2.0, and 
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both the imaginary and real parts of p2 vanish for a range of R values. Because the rest 
of the Floquet multipliers for e = 2.0 vanish (not shown here), all the dominant Floquet 
multipliers of e = 2.0 are related to the period-doubling bifurcation points. 

Figure shows the entrainability W and the period T for e = 2.5. We also see the rapid 
increase in W around R ~ 1.9, and the period T increases up to T = 1.3 for 1 < i? < 1.9 
and then decreases for R > 1.9. Figures [T2E and F show the first pi and the second p2 
dominant Floquet multipliers for e = 2.5, respectively, pi is a pure real multiplier with a 
minimum around R ~ 2.2, which does not coincide with the position of the maximum of W. 
Although the magnitude of p2 is smaller than pi, it gains a peak at i? ~ 1.9 identical to the 
peak position of W. Thus, the system is better entrained in the vicinity of the saddle-node 
bifurcation, and there exists an optimal coupling strength for the entrainability. 

We next show the entrainability W, the period T, and the Floquet multipliers Pi as a 
function of e for the relaxation oscillator. Figure [T3K is a dual-axis plot of the entrainability 
W and the period T ioi R = 1.3. We see very little dependence of the entrainability or 
the period on e. Figures [T3B and C show the first pi and the second p2 dominant Floquet 
multipliers, respectively; pi approaches 1 with decreasing e, indicating that the oscillator 
approaches the saddle-node bifurcation points, whereas p2 vanishes for a range of values of 
e. Figure [T3D shows the entrainability W and the period T as a function of e for R = 1.9, 
where the period T does not strongly depend on the coupling strength e. In contrast, the 
entrainability W exhibits nonmonotonic behavior and achieves the local maximum around 
e = 2.5. The first dominant Floquet multiplier pi in Figure fT3E shows that the oscillator 
approaches the period-doubling bifurcation when approaching e = 2.3 from above. Although 
the first dominant Floquet multiplier shows the period-doubling bifurcation, its magnitude 
does not qualitatively explain the existence of the entrainability peak around e = 2.5. In 
contrast, the second dominant Floquet multiplier p2 in Figure [T3F shows a positive-valued 
peak around e = 2.5, and the position agrees with that of the entrainability. In Figure [T31D. 
enhancement is only observed inside a narrow region where e < 2.8. 

4. DISCUSSION 

Mismatched periods enhance entrainability for both smooth and relaxation oscillators, 
and the enhancement is related to the dominant Floquet multipliers. The behavior as a 
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function of the model parameters {R and e) can be explained qualitatively by a pure real 
Floquet multiplier. It has been previously shown that in the vicinity of the bifurcation 
points, the oscillator is more sensitive to an external signal 47|, |48| . Entrainability is thus 
enhanced because the period- mismatch drives the coupled oscillator to the vicinity of saddle- 
node bifurcation. The pure real multiplier has a strong impact on the entrainment property 
(Figures [TOl [TT| [T2D-F and[T3D-F). as do the Floquet multipliers related to other bifurcation 
types (Figures [T2K-C). In Figure [T2]A. we see that even when neither the first nor the second 
dominant Floquet multiplier is a positive real multiplier, enhancement is induced by the 
period-doubling bifurcation points. 

From our analysis above, we saw that bifurcation is less likely to occur in a smooth 
oscillator than in a relaxation oscillator. Nevertheless, the enhancement of entrainability 
(an acute peak for a relaxation oscillator and an obtuse peak for a smooth oscillator) is seen 
in a wider region for a smooth oscillator than for a relaxation oscillator (Figures [TOK and 
[T2] D for R dependence, and Figures [TTlA and [T3t) for e dependence). Although a smooth 
oscillator may be better entrained. Figures [Eland [9] show that the amplitude of the oscillation 
is more strongly affected in the smooth oscillator. When increasing the period mismatch R, 
the amplitude of yi increases at the outset and then starts to decrease; the decrease starts 
sooner for the smooth oscillator than it does for the relaxation oscillator (see Figure H] and [6]) . 
Since limit-cycle oscillations with smaller amplitude are more sensitive to external stimuli, 
the reduction of the amplitude that is induced by the period mismatch is also a cause of the 
improved entrainability. However, the main contributing factor to the enhancement is the 
bifurcation on the limit cycle. This is because the behavior of the entrainability W agreed 
with that of the Floquet multiplier, which does not depend on the scaling of the values. 

As mentioned in the introduction, many circadian oscillators consist of multioscillatory 
networks, each of which has different components. As an example, early works on Gonyaulax 
reported periodic bioluminescence (24 h to 27 h: average 25 h) and aggregation (15 h to 
23 h: average 20 h) patterns {2, 3|. These two oscillators are synchronized in the presence 
of white light, whereas they are decoupled in the presence of dim red light. The ratio of 
the two periods is about 1.25 (see Figure 2(c) in Roenneberg j2|), which coincides with our 
results of i? ~ 1.4 for enhancement in the smooth oscillator (see Figure [TOK). 

We have shown that the oscillators are better synchronized to an external signal when they 
are close to bifurcation points (mainly at the saddle-node but also partially at the period- 
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doubling bifurcations). Regarding the period-doubling bifurcation in realistic systems, we 
note that experimental results showed the existence of a circabidian rhythm, with a 48 h 
period (double the circadian period). It has been reported that such circabidian rhythms 
can occasionally emerge under controlled conditions [49]. The bidian rhythm may result 
from a period-doubling bifurcation of a circadian rhythm that is close to a bifurcation point. 

Regarding stochasticity, gene expression is subject to two type of noise (stochasticity): 
one from the discrete molecular species (intrinsic noise) and the other from environmental 



variability (extrinsic noise) 50N53|. We may take into account stochastic effects in order to 



make our findings biochemically realistic. For example, a collection of identical oscillators 
may exhibit damped oscillation as the averaging effect of stochasticity Nevertheless, 
we employed a deterministic model (i.e., ordinary differential equations) in our analysis for 
two reasons. First, detailed stochastic modeling, such as a continuous-time discrete Markov 
chain, is not tractable: the only solution is an exhaustive, time-consuming calculation using 
the Monte Carlo method without any theoretical guarantees. Second, without analytical 
calculations, the theoretical arguments for the enhancement would be difficult. The en- 
hancement we detected at the bifurcation on the limit cycle (the saddle node bifurcation) is 
unlikely to be found by numerical simulations, but can be deduced from Floquet multipli- 
ers. Although our present analyses do not incorporate stochastic aspects, we consider that 
our results help in the understanding of the stochastic dynamics in modeling mismatched 
periods. We will consider such stochasticity in a future study. 

Another topic yet to be investigated is when three or more oscillators have mismatched 
periods. When considering more than two oscillators, besides the oscillatory model, the 
topology of the coupling is important. We can make an educated guess at the consequences 
of a period mismatch in such a multioscillator model: if period inconsistency drives the 
coupled system to the vicinity of bifurcation on the limit cycle, the period mismatch will 
enhances entrainability. As long as the overall connection of the multiple oscillators can 
be simplified as "two mismatched oscillators", the entrainability of such a component is 
expected to be enhanced. 

In summary, we have used phase reduction and Floquet multipliers to show that entrain- 
ability is enhanced by mismatched periods in both smooth and relaxation oscillators. We 
focused on two identical oscillators that were symmetrically coupled to each other by using a 
deterministic approach. In order to make our results more useful for real biological systems. 
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however, it would be important to consider situations of three or more oscillators, asymmet- 
ric coupling, or structurally heterogeneous cases with stochasticity. These extensions are 
left to our future studies. 
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